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On the basis of the lattice Boltzmann method we have done a numerical experiment of a forced 
turbulence in real space and time. Our new findings are summarized into two points. First in 
the analysis of the mean-field behavior of the velocity field using the exit-time statistics we have 
verified Kolmogorov's scaling and Taylor's hypothesis for the first time in the simulation for the 
Navier-Stokes equation. Second in the analysis of the intermittent velocity fluctuations using a 
non-equilibrium probability distribution function and the wavelet denoising we have clarified that 
' the coherent vortices sustain the power-law velocity correlation in the non-equilibrium state. 
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I. INTRODUCTION 



Understanding turbulence has not been attained over millenniums. However, new approaches to turbulence [1] have 
been invented successively. 

Recently turbulence data in nature and laboratories are analyzed from new view points. One [2] is based on the 
exit-time statistics and discusses Kolmogorov's scaling [1]. The other [3] on the basis of the Renyi-Tsallis statistics 
[4,5] discusses the probability distribution functions (PDF). 

Our study starts from the application of the above two schemes to our numerical simulation data and aims to obtain 
some new insight of the turbulence. 
I | Before reporting our study of the turbulence we briefly review here our strategy approaching intermittency problems 
' O ■ in other non-equilibrium systems. In the numerical study of the critical spin-state of the Ising model [6] we have 
observed that the degree of the non-Gaussianity of the Renyi-Tsallis PDF represents that of the non-equilibrium. On 
the other hand, in the numerical study of a supercooled liquid [7] near the glass transition we have observed that 
the spatial distribution of the coherently moving regions is well analyzed by the singularity spectrum of multifractal. 
Since the Renyi-Tsallis statistics can describe the systems with multifractality, the above two observations are closely 
related. 

In these two systems the correlation length of the fluctuations becomes divergently large and we can expect the 
scale-invariance. A scale-invariant system has multifractal nature. On the other hand, we can also expect the scale- 
invariance in the turbulent systems at high Reynolds number [1]. Thus we adopt the same strategy to intermittency 
. problem in turbulence. 

The next section describes our simulation method. In the following sections the simulation data are analyzed at 
the mean- field level at first and subsequently taking fluctuations into account. In the final section a brief summary is 
■ given, 
-i— » ' 

Sh ; II- LATTICE BOLTZMANN SIMULATION 

We simulate a forced turbulence in real space and time and adopt the lattice Boltzmann method [8,9] as one of the 
easiest ways for such a purpose. 

The distribution function fi(r,t) of the particle with the velocity Cj at position f and time t obeys the lattice 
Boltzmann equation [8,9] 



fi(f+ Ci, t + 1) - Mr, t) = -u[fi(f, t) - fl{r, t)], (1) 

where we have adopted the Bhatnagar-Gross-Krook (BGK) model for the collision term and the time step has been 
chosen as unity. The position r of the particle is restricted on the cubic lattice and the lattice spacing is chosen as 
unity. The local equilibrium distribution ff(r,t) is assumed to be given as 

/P(f, f) = w iP [l + 3@ • u) + |($ ■ uf - • u)]. (2) 

The density p and the velocity u of the fluid are given by the sum of the contributions of the particles as p — J2i fi{^ t) 
and pu = J2 i fi(f,t)ci, respectively. Using the multiscale analysis the Navier-Stokes equation, which is supposed to 
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be able to describe turbulence, is derived from the present BGK lattice Boltzmann equation and the corresponding 
viscosity v = (1/u) — l/2)/3. 

In the following we use a 15- velocity model and Ci is chosen as c\ = (0, 0, 0), &i = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) for 
i = 2, 3, ■ • •, 7 and c, = (±1, ±1, ±1) for i = 8, 9, — , 15. Then the weight factor Wi for the local equilibrium distribution 
is determined as w\ = 2/9, Wi = 1/9 for i = 2, 3, • • •, 7 and Wi — 1/72 for i = 8, 9, • • •, 15. The relaxation frequency lo 
is chosen as w — 1.94. 

In order to simulate forced turbulence we add the forcing term gi ■ F(f, t) due to an applied body force F(f, t) to 
the right hand side of Eq. (1). In order to reproduce the Navier-Stokes equation we should choose as $ = Cj/10 for 
the present 15-velocity model [8]. For simplicity we apply a solenoidal force in y-direction [9], F = (0, F y (z,t),0), 
and F v (z,t) is a Gaussian white noise whose variance o\f is a function of the lattice coordinate in z-dircction, 
<jp = |0.01 x sin(27rfc/z/L)|, where L is the linear dimension of the simulation region of cubic box. We have chosen 
as L = 200 so that 1 < x, y, z < 200. We have adopted the periodic boundary condition. Although we have chosen as 
kf = 4, the results of the simulation is insensitive to the choice. 

In the initial state of the simulation the fluid density is uniform, p = 1, and the spatial distribution of the fluid 
velocity u is chosen to be random. 

In Fig. 1 we show a snapshot of the vorticity field, u> = V x u, in order to visualize a turbulent structure. We see 
an inhomogcncous distribution of the vortices. 

The Reynolds number Re for our simulation [8,9] is estimated as Re = (u 2 )L/2tti' and time-dependent as shown 
in Fig. 2. Here (• • •} represents the spatial average. In our simulation we can realize a turbulent state for relatively 
small Reynolds number, since we add random force. 



III. MEAN-FIELD BEHAVIOR 



In this section we discuss the mean-field aspects of the velocity field. 

Turbulence is one of the typical phenomena with multiscalc motions. Each phenomena of a scale strongly couples 
with all the other scales of turbulent motion. In order to analyze such a system a scale-dependent entropy, the so-called 
e-entropy, works well [2]. For example, the time series of the velocity fluctuation in turbulence leads to a non-trivial 
scaling relation of the e-entropy, h(e) oc e~ 3 , expected from Kolmogorov's scaling [2]. The existing experimental data 
are consistent with this scaling [2] . In this section we try to show the consistency of the scaling with our numerical 
experiment simulating the Navier-Stokes equation. 

We focus our attention to the time series of the y-componcnt u y (f,t) of the fluid velocity u(r,t) as shown in Fig. 
3, since our system is anisotropic due to the forcing. This turbulent signal is characterized by the exit-time r(e; r) at 
position f. At the exit-time t = r(e;r) the first exit satisfying the condition, \u y (r, to + t) — u y (f,to)\ > e/2, occurs. 
Here the time r is measured from to and to = WOO in our simulation regarding the states at the first 1000 time-steps 
as transitional. 

The mean of the exit-time, r(e) = (r(e; r)), is related to the e-entropy, h(e), as [2] 

h(e) oc -L (3) 

In the region where Kolmogorov's scaling, the 4/5 law [1,2], holds the velocity difference between two points behaves 
as 

(\u(r + R,t )-u(f,to)\)^\R\ 1/3 . (4) 

Using Taylor's hypothesis [1,2], we obtain 

(\u(f,t + T) -u(r,t )\) cxr 1 / 3 , (5) 

as a mean-field description. In spite of the absence of the global flow in our simulation Taylor's hypothesis is valid in 
a local sense [11] as shown in the following. Namely in the mean-field description the velocity difference in temporal 
and spatial directions have the same fractal scaling exponent. With Eq. (5) the definition of the exit-time leads to 
r(e) 1 / 3 oc e. From Eq. (3) we can conclude that h(e) oc e 3 in the above mean-field description [2]. 

In accordance with the above scaling analysis our simulation data [10] for l/r(e) shown in Fig. 4 have the scaling 
region for e > e where e ~ 0.04. Thus by our simulation we have established that h(e) oc e~ 3 . For e <C e the exit 
in the simulation occurs within one time step and the curve saturates at small e. These behaviors are consistent with 
the experimantal data [2]. 
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In Fig. 5 the result of the similar analysis in spatial direction is shown. In place of r(e) we measure r(e) where at 
the exit-length r — r(e; r) the first exit satisfying the condition, \u y {x + r, y, z, t) — u y (x, y, z,t)\ > e/2, occurs at a 
time t. Here r(e) = (r(e;r)). Because of the smallness of the number of the available data in space direction, 200, 
in comparison with that in time direction, 7000, the scaling region is narrow. Since l/r(e)and l/r(e) have the same 
scaling exponent, Taylor's hypothesis is confirmed to be valid. 

In conclusion we have confirmed Kolmogorov's scaling and Taylor's hypothesis at the same time. In contrast to 
the evaluation of the energy spectrum, which is usually employed for testing Kolmogorov's scaling and determined by 
the two-point correlation function, Kolmogorov's scaling is easily observed in the e-entropy, since it is a mean-field 
description and fluctuations are averaged out. Such a mean-field scaling is a unifractal description and fluctuations 
can be taken into account in a multifractal analysis as shown in the following. 



We have discussed the mean-field aspects of the velocity field in the preceding section. Next we discuss the 
fluctuations of the velocity field. In this section we make a qualitative observation using the correlation maps. 

In Fig. 6 we show the space correlation. The correlation decays as the length scale is increased. In Fig. 7 we show 
the time correlation. The correlation decays as the time scale is increased. At small length or time scale the strong 
correlation leads to a non-equilibrium state which is identified as the coherent vortex in the following section. 

We see a multiscale structure in the correlation map, Fig. 8, which is one of the characteristics of turbulence neither 
periodic nor random. 



We have discussed the qualitative aspects of the velocity fluctuation in the preceding section. Next we try to 
quantify the fluctuation by two measures. One is the non-equilibrium parameter used in this section and the other is 
the singularity spectrum in the next section. 

In order to analyze the distribution of the fluctuation we use the Renyi- Tsallis PDF [4,5] 



for a variable X. Here the parameter q represents the degree of non-equilibrium [6]. To be precise, the deviation of q 
from unity is the measure for the degree of non-equilibrium, since the equilibrium Gaussian distribution is expressed 
by Eq. (6) in the limit of q — > 1. 

The PDF in Eq. (6) is derived from the Renyi or Tsallis entropy by using the maximum entropy principle [12]. 
The Tsallis entropy is non-extensive, while the Renyi entropy is extensive. The functional form of the distribution is 
independent of the choice of the entropy, the Renyi or Tsallis, so that our result is independent of the extensivity of 
the entropy. 

We show the PDF for the velocity difference between two points in Fig. 9 and the vorticity difference in Fig. 10. 
The q- value of the PDF is length-scale dependent. Similarly the PDF for the velocity or vorticity difference between 
two times at a point is time-scale dependent. In any case the degree of non-equilibrium is larger at smaller spatio- 
temporal scales where the correlation survives. As shown later on, the formation of the coherent vortex structure 
sustains the strong correlation leading to non-Gaussian PDF. At larger scales the contribution of the structureless 
incoherent background among vortices dominates. This incoherency leads to Gaussian PDF. The Gaussian PDF is a 
parabola in the semi-logarithmic plane as Fig. 9 and Fig. 10. 

Although our data in Fig. 9 are not enough to discuss the scaling between q and r quantitatively, we have to derive 
such a relation in our future study. Already a scaling relation between q and r has been proposed [13] under the 
assumption of a cascade picture. On the other hand, another scheme with a r-independent g-value has been proposed 
[12] under the assumption of another cascade picture. In contrast to these two approaches our analysis has nothing 
to do with cascades. 

In Fig. 11 we show the PDF for the vorticity itself taken from a snapshot of the vorticity field. This PDF is 
also non-Gaussian. The Gaussian components can be removed from the snapshot by using the wavelet denoising 
technique [14]. It can be seen from the coincidence between the left and right panels in Fig. 12 that the non-Gaussian 
components corresponds to the vortices. Here the Gaussian component is filtered out by the wavelet denoising in 
the right panel. This visualization is not new, since it has already been done in the study of the coherent vortex 
simulation [14]. Our new contribution is to clarify the origin of the scale-dependent q- value in terms of the coherent 



IV. FLUCTUATIONS: QUALITATIVE OBSERVATION 



V. FLUCTUATIONS: NON-EQUILIBRIUM NATURE 
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vortex structure. The strong correlation leading to the non-Gaussian PDF at smaller scales is sustained in the inner 
region of the vortices, while the PDF at larger scales than the vortex size tends to Gaussian. In other words the inner 
region is coherent and the outer is incoherent. 

VI. FLUCTUATIONS: MULTIFRACTAL NATURE 

We have discussed the local nature of the vortex as the elemntary excitation in the preceding section. In this section 
we discuss the spatial distribution of the vortices in a global point of view. 

First we make a qualitative observation. Already as seen in Fig. 12 the spatial distribution of the coherent vortices 
is intermittent. In Fig. 13 we also see a intermittent behavior for the derivative of the vorticity. Similar intermittent 
behavior is observed in experimental data [15]. 

Next we quantify the intermittency by introducing the measure, ^k(l), defined as the probability of finding a 
coherent vortex in the fc-th cubic box of the linear dimension I in a snapshot of the vorticity field. The total number 
of the boxes Nl = (L/l) x (L/l) where L is the linear dimension of the simulation region. The measure behaves as 

Wb(0«J a *. (7) 
for small l/L. The number density N(a) for the exponent a k is defined as 

N(a) =J2 6 (a- a k ) ocl- f ^\ (8) 
k 

This relation defines the singularity spectrum, f(a), which is the key quantity to qualify multifractal systems. Using 
the normalized q-th moment fj,k(q,l) of the measure /Zfc(Z), 

N L 

wfe(«,0 = {/*(0}7£W(0}*, (9) 
k'=i 

the sigularity spectrum f(a) is obtained by the following formulae, [16] 

N L 

a(q) - J2 InMO/ H^/L), (10) 
k=i 

and 

N L 

/(«(«)) = E M<z> ln Ms> 0/ Hi/L). (n) 

*;=i 

In Fig. 14 we show the singularity spectrum for the coherent vortex. For simplicity we have used the definition 
Nk/N^ = l ak where N k is the number of lattice points satisfying the condition io z > (u> z ). The spectrum is similar to 
the typical one measured for the energy-dissipation rate by experiments [15] or simulations [17]. 

In the last part of this section we give some speculations. As seen in our previous study [7] the width of the 
singularity spectrum depends on the degree of intermittency so that we expect some Reynolds-number dependence 
of the spectrum. It has been discussed by many authors [18] and should be claryfied in future systematic study. 
In our present study the density of the vortex is dilute as seen in Fig. 12. Thus the singularity spectrum in Fig. 
14 describes the spatial distribution of relatively free vortices. As the Reynolds number is increased the density 
increases [19]. In this case the interaction among vortices becomes important and the correlation length of the 
vorticity fluctuation becomes large. While we can observe the vortex only as an individual elementary excitation 
in our numerical experiment, some collective excitation is expected to dominate at higher Reynolds number. In the 
limit of divergently large Reynolds number the correlation length becomes divergently large so that we can expect full 
scale-invariance. We can find a resemblance to the case of the scaling theory [20] in polymers where an ideal scaling 
relation is realized for dense solutions where polymers are strongly entangled. In the limit of high Reynolds number 
each boxes counting the coherent vortex will be filled by almost equal number of vortices so that intermittency will 
disappear and unifractal Kolmogorov's scaling will prevail. 



4 



VII. SUMMARY 



Numerical simulation data, in real space and time, for a forced turbulence on the basis of the lattice Boltzmann 
method have been analyzed by unifractal and multifractal schemes. 

Our new findings are summarized into two points. First in the unifractal analysis using the exit-time statistics we 
have verified Kolmogorov's scaling and Taylor's hypothesis at the same time. Second in the analysis using the Renyi- 
Tsallis PDF and the wavelet denoising we have clarified that the coherent vortices sustain the power-law velocity 
correlation in the non-equilibrium state. 

Finally in the multifractal analysis it is clarified that the intermittent distribution of the coherent vortices in 
space-time is described as a multifractal. 

This work was supported in part by a Grand-in-Aid for Scientific Research from the Ministry of Education, Culture, 
Sports, Science and Technology of Japan. 
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FIG. 1. The snapshot of the ^-component of the vorticity, u) Z) in the £y-plane at z = 100 and t — 4000. 

FIG. 2. The Reynolds number Re as a function of time t. 

FIG. 3. The {/-component of the fluid velocity u y as a function of time t at f = (100, 100, 100). 

FIG. 4. The inverse of the mean exit-time r(e). The exits are measured at 200 a lattice points for 1000 <t< 8000. Here the 
dots and straight line represent the simulation data and e~ 3 law, respectively. 

FIG. 5. The inverse of the mean exit-length r(e). The exits are measured at 200 lattice points at t = 6000. Here the dots 
and straight line represent the simulation data and e -3 law, respectively. 

FIG. 6. The space correlation of the velocity fluctuation at t = 4000. The horizontal and vertical axes represent u y (x, y, z, t) 
and u y (x,y + y C orr, z,t), respectively. The space difference y COII is taken as y COII — 1 for the left, y corr — 4 for the center and 
J/corr = 30 for the right. 
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FIG. 7. The time correlation of the velocity fluctuation at r = (100, 100, 100). The horizontal and vertical axes represent 
u v {r,t) and u y (f,t + t COII ), respectively for 1000 < t < 8000. The time difference t COII is taken as t COII — 10 for the left, 
icorr = 100 for the center and t COTT — 500 for the right. 

FIG. 8. The correlation map for the velocity difference between two times ti and tj for 1000 < to + ti < 8000 and 
1000 < t Q + tj < 8000 with t = 1000. Namely < i < 7000 (horizontal axis) and < j < 7000 (vertical axis). The 
dot is marked at if D(i,j) < 8 where D(i,j) = \u y (f,ti) — u y (f,tj)\ and 9 = 0.2D max with D m „ being the maximum of 

D(i, j). The position f is fixed at (100, 100,50). 

FIG. 9. The unnormalized PDF P(X) for the velocity difference X = \u y (x, y + Ay, z,t) — u y (x,y, z,t)\ between two points 
(x,y + Ay,z) and (x,y,z) at t = 4000. The difference in y-direction is chosen as Ay — 1,4,30. The data are fitted by the 
Renyi-Tsallis PDF with q = 1.218, 1.164, 1.100 for Ay = 1,4, 30, respectively. 

FIG. 10. The unnormalized PDF P(X) for the vorticity difference X = \io y (x, y + Ay, z, i) — u y (x, y, z, t)\ between two points 
(x,y + Ay,z) and (x,y,z) at t — 6000. The difference in y-direction is chosen as Ay — 1,4,30. The data are fitted by the 
Renyi-Tsallis PDF with q = 1.240, 1.209, 1.179 for Ay = 1, 4, 30, respectively. 

FIG. 11. The unnormalized PDF P(X) for the vorticity X = \u y (x, y, z, t)\ at t = 6000. The data are fitted by the 
Renyi-Tsallis PDF with q = 1.392. 

FIG. 12. The snapshot of the z-component of the vorticity, u) z , in the xy-pl&ne at z — 100 and t = 4000. In contrast to Fig. 
1 the regions satisfying the condition, cu z > (u) z ) + 2a = O, are depicted in the left panel. Here (u) z ) is the spatial average and 
a is the standard deviation. In the right panel the high vorticity regions satisfying the same condition, Cu z > O, are depicted 
where u> z is the non-Gaussian component after the wavelet denoising filtering out the Gaussian component. By the restriction 
of the algorithm of the wavelet transform we have used only 128x128 lattice points. 

FIG. 13. The spatial distribution of the vorticity gradient, Auj z = u z (x + 1, y, z, t) —cj z (x, y, z, t), in a;-direction (1 < x < 200) 
with y = 2 = 100 at t = 6000. 

FIG. 14. The singularity spectrum of the coherent vortex with l/L = 1/10. The data correspond to —90 < q < 95. 
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